module lcmfd2g !larsen's CMFD
    use mtbl2g
!    use geom, only: nxy, nz
    implicit none
    real(8), private, pointer, dimension(:,:,:)     :: phi
    real(8), private, pointer, dimension(:,:)       :: psi, psip
    type(mb2g), private :: mtcmfd, mtcmfdorg
    integer, private :: maxout, maxin 
    real(8), private :: ccout, ccin
    integer, private :: nx, ny, nxy, nz
    logical, private :: l1gcmfd, ifwldt
    real(8), private, pointer, dimension (:) :: h2x, h2y
    real(8), private, pointer, dimension (:,:,:) :: phi1Lsx, phi1Lsy
    real(8), private, pointer, dimension (:,:,:) :: trlx, trly
    real(8), private, pointer, dimension(:,:,:,:)   :: dtilde, dhat
    real(8), private :: k_eff
    integer, private :: SGDHAT(4) ! sign of dhat
    data SGDHAT / -1, 1, 1, -1/
    logical, private, pointer, dimension(:,:,:,:)       :: vacmesh
    real(8), private, pointer, dimension(:,:)           :: psiacc
	integer, private, pointer, dimension(:)				:: sx, sy, ex, ey   
contains

subroutine resetpsi_lcmfd2g
    psiacc=0;
endsubroutine

subroutine initlcmfd2g
    use allocs
    use bicg2g, only : mallocbicg2g
    implicit none
    include "pspec.FI"
    include "xsec.FI"
    integer :: i,j

    nx=ncmx;ny=ncmy;nxy=nx*ny;nz=1
    
    call dmalloc0(h2x, 0, nx)
    call dmalloc0(h2y, 0, ny)
    
    h2x(0)=hcmx/2
    h2x(nx)=hcmx/2
    do i=1, nx-1
        h2x(i)=hcmx
    enddo 

    h2y(0)=hcmy/2
    h2y(ny)=hcmy/2
    do i=1, ny-1
        h2y(i)=hcmy
    enddo 
    
    call dmalloc0(phi1lsx,  1,ng2, 0,nx, 0,ny)
    call dmalloc0(phi1lsy, 1,ng2, 0,nx, 0,ny)
    call dmalloc(trlx, ng2, nx, ny)
    call dmalloc(trly, ng2, nx, ny)
    call dmalloc0(dtilde, 1,2, 1,ng2, 0,nx, 0,ny)
    call dmalloc0(dhat,   1,2, 1,ng2, 0,nx, 0,ny)
    call dmalloc(vacmesh, ng2, nx, ny, 1)
    
    
!    nx=ncmx;ny=ncmy;nxy=nx*ny;nz=1
    call dmalloc0(phi,1,2,0,nxy,0,1)
    call dmalloc(psi,nxy,1)
    call dmalloc(psip,nxy,1)
    call newmb2g(mtcmfd,nxy,nx,ny, 1)
    call newmb2g(mtcmfdorg,nxy,nx,ny, 1)
    call dmalloc(psiacc, nxy, 1)

!    call mallocbicg2g(nxy,nx,ny,nz)
    ifwldt = lwldt
!    if (ifwldt) print *, "Wielandt acceleration turned on!"
    if (gc2(ng) .eq. 1) l1gcmfd = .true.
    
    psip=1;

    maxout = max_out
    maxin  = max_in
    ccout = err_out
    ccin  = err_in
!    phi=1.

	call dmalloc(sx,ny)
	call dmalloc(ex,ny)
	call dmalloc(sy,nx)
	call dmalloc(ey,nx)
	do j=1, ny
		do i=1, nx
			if (sx(j) .eq. 0 .and. map_exp2(i,j) .ne. -1) sx(j)=i
			if (map_exp2(i,j) .ne. -1 ) ex(j)=i
		enddo
	enddo
	do i=1, nx
		do j=1, ny
			if (sy(i) .eq. 0 .and. map_exp2(i,j) .ne. -1) sy(i)=j
			if (map_exp2(i,j) .ne. -1 ) ey(i)=j
		enddo
	enddo
!	sx=1;ex=nx;sy=1;ey=ny
endsubroutine

subroutine sollcmfd2g
    integer :: iter
    real(8) :: k, err, dr, i, g
    logical :: iflog
!    ifwldt = .false.
    iflog  = .true.
    k=1._8
    do i=1, nxy
        do g=1, ng2
            phi(g,i,1) = 1._8 !drand(0)
        end do
    end do
    call caldtilde
    call makemt(.false.)
    call poweriter(k, iter, err, dr, .false., iflog)
!    call poweriter(k, iter, err, dr, .true.)
endsubroutine

subroutine lcmfd2gacc(k_in, err, iter, psidif, onfdb, powout)
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    include "statistics.FI"
!    ! for debugging ========================
!    include "hmgz.FI"
!    ! ======================================
    integer :: iter
    real(8), intent(inout) :: k_in
    real(8) :: k, err, dr, g, sumadj
    real(8), intent(out) :: psidif
    real(8) :: upper, lower
    logical :: ifdhat, iflog
    integer :: i,j,l
    integer :: nfispin, nzero
    real(8) :: nrmf, tpsi
    logical, intent(in) :: onfdb
    real(8), dimension(:,:) :: powout
    
    ifdhat = .true.
    iflog  = .false.
    k=k_in*1.1
    phi=1.
    
	do g=1, ng2
        l=0;
        do j=1, ny
            do i=1, nx
                l=l+1
                if (map_exp2(i,j) .eq. -1) cycle
                phi(g,l,1)=hphi(g,i,j)
                if (hphi(g,i,j) == 0._8) then
                    htot(g,i,j)=htot2(g,i,j)
                    hdif(g,i,j)=1._8/(3*htot2(g,i,j))
                    hnuf(g,i,j)=hnuf2(g,i,j)
                    hrmv(g,i,j)=hRmv2(g,i,j)
                    hmsct(:,g,i,j)=hmsct2(:,g,i,j)
                    vacmesh(g,i,j,1)=.true.
                    nzero=nzero+1
                else
                    vacmesh(g,i,j,1)=.false.
                endif
            enddo
        enddo
    enddo
    call caldtilde
    call caldhat
    if (l1gcmfd) then ! the one-group CMFD
        call makemt1g(ifdhat)
    else
        call makemt(ifdhat)
    endif
    call poweriter(k, iter, err, dr, ifwldt, iflog)
    if (ifwldt .eq. .false.) print *, dr
    ! check the solution
    nfispin=0
    do l=1, nxy
        if ( psi(l,1)<0) then
            print *, "negative flux"
			k=0.0_8
            iter = 0
            err  = 1.0_8
            do i=1, nxy
                phi(:,i, 1) = 1._8
            enddo
            fdbsafe = .false.
            return
        elseif (psi(l,1)>0) then
            nfispin = nfispin+1
        endif
    enddo
    
    
    
    ! error check 
    nrmf =1._8/sum(psi)
    psi=psi*nrmf

    psiacc=psiacc+psi
    nrmf = 1._8/sum(psiacc)
    ! print 
!    open (123, file = "cmfdpsi.txt", status='unknown')
!    do j=ny, 1, -1
!        write(123, '(2000(es12.5,1x))'), psi((j-1)*nx+1:j*nx,1)
!    enddo
!    close(123)
    
    upper=0;lower=0;
    if (ireft .eq. 1) then
        if (hybridcmfd) then
            do l=1, nxy
                upper = upper + (refcmfd(l)-psiacc(l,1)*nrmf)**2
            enddo
        else
            do l=1, nxy
                upper = upper + (refcmfd(l)-psi(l,1))**2
    !            lower = lower + refcmfd(l)
            enddo
        endif
        lower=1
    else
    !    nfispin=0;
        do l=1, nxy
            upper = upper + (psi(l,1)-psip(l,1))**2
            lower = lower + psi(l,1)
    !        nfispin = nfispin+1
        enddo
        psip=psi
    endif
    do j=1, ny
        powout(:,j)=psi((j-1)*ny+1:j*ny,1)
    enddo
    
    psidif=sqrt(upper*npowasmref)/lower !*sqrt(1.0_8*nfispin)
    fdbsafe = .true.

    nrmf=nrmf*nfispin
    psi = psi * nrmf *vcm
    phi = phi * nrmf
     
    k_in = k
    k_eff = k
    l = 0
    do j=1,ny
    do i=1, nx
        l = l+1
        if (map_exp2(i,j) .eq. -1) then ! for vacuum assembly
            adjp(i,j) = 0._8
        else
            adjp(i,j)=psi(l,1)
            sumadj = sumadj + adjp(i,j)
!           fission weight map
            if (psi(l,1)>0) fiswtmap(i,j)=1._8/psi(l,1)
        endif
    enddo
    enddo
    
!    adjp =1
!    fiswtmap =1
    
    do j=1,ny
    do i=1, nx
        adjp(i,j) = adjp(i,j)/sumadj
    enddo
    enddo

    if (unifsrc .and. onfdb) then
        psi=psi*nfispin/sum(psi)
        l=0
        do j=1, ncmy;do i=1, ncmx
            l=l+1
            pwadj(i,j) =1._8/psi(l,1)
            pwadjr(i,j)=psi(l,1)
        enddo;enddo
    endif

!    if (hybridcmfd ) then
!        do i=1, nx*ny
!            tpsi=psi(l,1)
!            stpsia(l,1)=stpsia(l,1)+tpsi
!            stpsiasq(l,1)=stpsiasq(l,1)+tpsi**2
!            l=l+1
!        enddo
!    endif

!    do j=1,ny
!        adjp(:,j)=psi((j-1)*nx+1: j*nx,1)
!    end do
!    adjp = adjp/sum(adjp)    
endsubroutine

subroutine caldtilde
    use define
    implicit none
    
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    
    integer i, j, k, g, oi, oj, bt
    real(8)    betaL, betaR, compL, compR
    
    Dtilde = 0.
    do g=1, ng2
        do j=1, ny
            i=sx(j)-1
            dtilde(1,g,i,j)=1._8/3/(htot(g,i+1,j))/h2x(i)
            do i=sx(j), ex(j)-1
                if (map_exp2(i,j) .eq. -1) cycle
                dtilde(1,g,i,j)=1._8/3*((hcmx+hcmx)/(htot(g,i,j)*hcmx+htot(g,i+1,j)*hcmx))/h2x(i)
            enddo
            dtilde(1,g,i,j)=1._8/3/(htot(g,i,j))/h2x(i)
        enddo
    enddo

    do g=1, ng2
        do i=1, nx
            j=sy(i)-1
            dtilde(2,g,i,j)=1._8/3/(htot(g,i,j+1))/h2y(j)
            do j=sy(i), ey(i)-1
                if (map_exp2(i,j) .eq. -1) cycle
                dtilde(2,g,i,j) = 1._8/3*((hcmy+hcmy)/(htot(g,i,j)*hcmy+htot(g,i,j+1)*hcmy))/h2y(j)
            enddo
            dtilde(2,g,i,j)=1._8/3/(htot(g,i,j))/h2y(j)
        enddo
    enddo
end subroutine

subroutine caldhat 
    use define
    implicit none
    include "xsec.FI"
    include "pspec.FI"
    include "arrays.FI"
    include "hmgz.FI"
    integer :: i, j, k, g, oi,oj, mi, mj
    real(8) :: bt, current
    real(8) :: hcm(4)
    real(8) :: check
    
    real(8) :: jsb, phi1sb, jst, js2t, l1st, dt, dt2
    
    
!    open(3001, file="dhatx.txt", status="unknown")
!    open(3002, file="dhaty.txt", status="unknown")
    hcm(1:2) = hcmx
    hcm(3:4) = hcmy
    
    do g=1, ng2
        do j=1, ny
            i=sx(j)-1
            dt=dtilde(1,g,i,j)
            dt2=dt*3/2
            phi1Lsx(g,0,j)=dt2*hcmx*(trcwx(g,i+1,j))
            do i=sx(j), ex(j)-1
                dt=dtilde(1,g,i,j)
                dt2=dt*3/2
                phi1Lsx(g,i,j)=dt2*hcmx*(trcwx(g,i,j) + trcwx(g,i+1,j))
            enddo
            dt=dtilde(1,g,i,j)
            dt2=dt*3/2
            phi1Lsx(g,i,j)=dt2*hcmx*(trcwx(g,i,j))
        enddo
    enddo

    do g=1, ng2
        do i=1, nx
            j=sy(i)-1
            dt=dtilde(2,g,i,j)
            dt2=dt*3/2
            phi1Lsy(g,i,0)=dt2*hcmy*trcwy(g,i,j+1)
            do j=sy(i), ey(i)-1
                dt=dtilde(2,g,i,j)
                dt2=dt*3/2
                phi1Lsy(g,i,j)=dt2*hcmy*(trcwy(g,i,j) + trcwy(g,i,j+1))
            enddo
            dt=dtilde(2,g,i,j)
            dt2=dt*3/2
            phi1Lsy(g,i,j)=dt2*hcmy*trcwy(g,i,j)
        enddo
    enddo
    
    do g=1, ng2
        do j=1, ny
            do i=1, nx
				if (map_exp2(i,j) .eq. -1) cycle
                trlx(g,i,j)=(phi2trx(g,i,j)-phi2trx(g,i,j-1))/hcmy*hcmx
                trly(g,i,j)=(phi2try(g,i,j)-phi2try(g,i-1,j))/hcmx*hcmy
            enddo
        enddo
    enddo
    
    
!    phi2x=0;phi2y=0;phi1lsx=0;phi1lsy=0
!    trlx=0;trly=0
    do g=1, ng2
        do j=1, ny
            i=sx(j)-1
            dt=dtilde(1,g,i,j)
            dt2=dt*3/2
            jsb    =-hj(1,g,1,j)
            phi1sb =phi1Lsx(g,i,j)
            jst    = -dt  *(hphi(g,i+1,j)-0)
            js2t   = -dt2 *(phi2x(g,i+1,j)-phi2x(g,i,j))
            l1st   =  dt2 *(trlx(g,i+1,j)+0)
!            phi1sb=0;js2t=0;l1st=0;
            dhat(1,g,i,j)=(jsb-phi1sb-(jst-js2t+l1st))/(hphi(g,i+1,j)+0)
!            write(3001, '(2000(e15.5, 1x))'), dhat(1,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
!            if (jsb .eq. 0._8) dhat(1,g,i,j)=0._8
            do i=sx(j), ex(j)-1
                dt=dtilde(1,g,i,j)
                dt2=dt*3/2
                jsb    =hj(2,g,i,j)
                phi1sb =phi1Lsx(g,i,j)
                jst    = -dt  *(hphi(g,i+1,j)-hphi(g,i,j))
                js2t   = -dt2 *(phi2x(g,i+1,j)-phi2x(g,i-1,j))
                l1st   =  dt2 *(trlx(g,i+1,j)+trlx(g,i,j))
!            phi1sb=0;js2t=0;l1st=0;
                dhat(1,g,i,j)=(jsb-phi1sb-(jst-js2t+l1st))/(hphi(g,i+1,j)+hphi(g,i,j))
!                if (jsb .eq. 0._8) dhat(1,g,i,j)=0._8
!                write(3001, '(2000(e15.5, 1x))'), dhat(1,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
            enddo
            dt=dtilde(1,g,i,j)
            dt2=dt*3/2
            jsb    =hj(2,g,i,j)
            phi1sb =phi1Lsx(g,i,j)
            jst    = -dt  *(0-hphi(g,i,j))
            js2t   = -dt2 *(phi2x(g,i,j)-phi2x(g,i-1,j))
            l1st   =  dt2 *(0+trlx(g,i,j))
!            phi1sb=0;js2t=0;l1st=0;
            dhat(1,g,i,j)=(jsb-phi1sb-(jst-js2t+l1st))/(0+hphi(g,i,j))
!            if (jsb .eq. 0._8) dhat(1,g,i,j)=0._8
!            write(3001, '(2000(e15.5, 1x))'), dhat(1,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
        enddo
    enddo
    
!    close(3001)
    do g=1, ng2
        do i=1, nx
            j=sy(i)-1
            dt=dtilde(2,g,i,j)
            dt2=dt*3/2
            jsb    =-hj(4,g,i,1)
            phi1sb =phi1Lsy(g,i,j)
            jst    = -dt  *(hphi(g,i,j+1)-0)
            js2t   = -dt2 *(phi2y(g,i,j+1)-phi2y(g,i,j))
            l1st   =  dt2 *(trly(g,i,j+1)+0)
!            phi1sb=0;js2t=0;l1st=0;
            dhat(2,g,i,j)=(jsb-phi1sb-(jst-js2t-l1st))/(hphi(g,i,j+1)+0)
!            if (jsb .eq. 0._8) dhat(2,g,i,j)=0._8
!            write(3002, '(2000(e15.5, 1x))'), dhat(2,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
            do j=sy(i), ey(i)-1
                dt=dtilde(2,g,i,j)
                dt2=dt*3/2
                jsb    =hj(3,g,i,j)
                phi1sb =phi1Lsy(g,i,j)
                jst    = -dt  *(hphi(g,i,j+1)-hphi(g,i,j))
                js2t   = -dt2 *(phi2y(g,i,j+1)-phi2y(g,i,j-1))
                l1st   =  dt2 *(trly(g,i,j+1)+trly(g,i,j))
!            phi1sb=0;js2t=0;l1st=0;
                dhat(2,g,i,j)=(jsb-phi1sb-(jst-js2t-l1st))/(hphi(g,i,j+1)+hphi(g,i,j))
!                if (jsb .eq. 0._8) dhat(2,g,i,j)=0._8
!                write(3002, '(2000(e15.5, 1x))'), dhat(2,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
            enddo
            dt=dtilde(2,g,i,j)
            dt2=dt*3/2
            jsb    =hj(3,g,i,j)
            phi1sb =phi1Lsy(g,i,j)
            jst    = -dt  *(0-hphi(g,i,j))
            js2t   = -dt2 *(phi2y(g,i,j)-phi2y(g,i,j-1))
            l1st   =  dt2 *(0+trly(g,i,j))
!            phi1sb=0;js2t=0;l1st=0;
            dhat(2,g,i,j)=(jsb-phi1sb-(jst-js2t-l1st))/(0+hphi(g,i,j))
!            write(3002, '(2000(e15.5, 1x))'), dhat(2,g,i,j), jsb-jst, jsb-phi1sb-(jst-js2t+l1st), jsb, phi1sb, jst, js2t, l1st
!            if (jsb .eq. 0._8) dhat(2,g,i,j)=0._8
        enddo
    enddo
    

end subroutine

subroutine makemt(ldh)
    use const, only :indm24
    use define
!    use geom
    implicit none    
    
    include "pspec.FI"
    include "hmgz.FI"
    
    logical, intent(in) :: ldh
    integer :: i, j, k, l, lo, g, gp, dr, d, nxs, nxe, si, sj, oi,oj
    real(8) ::offsum, val, dt, dh
    real(8) :: hcm(4)
   
    hcm(1:2) = hcmx
    hcm(3:4) = hcmy
    do k=1, 1
    l=0
    do j=1, ny
        nxs = 0
        nxe = 1
        do i=1, nx
            l=l+1
            if (map_exp2(i,j) .eq. -1) cycle ! for vacuum assembly
            ! check the first and the last mesh in a row
            if (nxs .eq. 0)  nxs = i
            nxe = i
            ! increment of real mesh index
            ! map for inner solver
            call setnodel(mtcmfd,i,j,l)
            !call setnxse(mtcmfd,j,1,nx)
            do g=1, ng2
                call setmb2g(mtcmfd,l,k,0,indm24(g,g),    hRmv(g,i,j)*vcm)
                do dr=1, 4
                    ! neighbor
                    lo = neighc(i,j,dr)

                    si=i;sj=j
                    if (dr>2) then
                        sj=j+(LTSIDXY(dr)-1)/2
                    else
                        si=i+(LTSIDXX(dr)-1)/2
                    endif
                    call neighc2(i,j,k,oi,oj)
                    if (oi>0 .and. oj>0) then
					    if (vacmesh(g, oi,oj,1) .or. vacmesh(g, i,j,1) ) dhat(LTSURFDIR(dr),g,si,sj)=0._8
					else
					    if (vacmesh(g, i,j,1) ) dhat(LTSURFDIR(dr),g,si,sj)=0._8
					endif
                    if (ldh) then ! with Dhat
                        dt=dtilde(LTSURFDIR(dr),g,si,sj)
                        dh=dhat(LTSURFDIR(dr),g,si,sj)
                        val=(-dt+dh*SGDHAT(dr))*hcm(dr)
                        if (lo .ne. 0) call setmb2g(mtcmfd,l,k,dr,g, val,lo)
                        val=(dt+dh*SGDHAT(dr))*hcm(dr)
                        call accmb2g(mtcmfd,l,k,0,indm24(g,g),val)
                    else    ! without Dhat
                        val=-dtilde(LTSURFDIR(dr),g,si,sj)*hcm(dr)
                        if (lo .ne. 0) call setmb2g(mtcmfd,l,k,dr,g, val,lo)
                        val=dtilde(LTSURFDIR(dr),g,si,sj)*hcm(dr)
                        call accmb2g(mtcmfd,l,k,0,indm24(g,g),val)
                    endif
                enddo
            enddo
            ! group coupling
            call setmb2g(mtcmfd,l,k,0,indm24(1,2),-hmSct(1, 2, i,j)*vcm)
            call setmb2g(mtcmfd,l,k,0,indm24(2,1),-hmSct(2, 1, i,j)*vcm)
        enddo ! loop for x-direction
        ! set the real meshes in a row
        call setnxse(mtcmfd,j,nxs,nxe) 
    enddo ! loop for j-directoin
    enddo ! loop for plane
    
    call copymb2g(mtcmfdorg, mtcmfd)
    call facilu2g(mtcmfd)
end subroutine

subroutine makemt1g(ldh)
    use const, only :indm24
    use define
!    use geom
    implicit none    
    
    include "pspec.FI"
    include "hmgz.FI"
    
    logical, intent(in) :: ldh
    integer :: i, j, k, l, lo, g, gp, dr, d, nxs, nxe, si, sj
    real(8) ::offsum, val
    real(8) :: hcm(4)
   
    hcm(1:2) = hcmx
    hcm(3:4) = hcmy
    do k=1, 1
    l=0
    do j=1, ny
        nxs = 0
        nxe = 1
        do i=1, nx
            ! increment of real mesh index
            l=l+1
            if (map_exp2(i,j) .eq. -1) cycle ! for vacuum assembly
            ! check the first and the last mesh in a row
            if (nxs .eq. 0)  nxs = i
            nxe = i
            ! map for inner solver
            call setnodel(mtcmfd,i,j,l)
            do g=1, ng2
                call setmb2g(mtcmfd,l,k,0,indm24(g,g),    hRmv(g,i,j)*vcm)
                do dr=1, 4
                    ! neighbor
                    lo = neighc(i,j,dr)
                    si=i;sj=j
                    if (dr>2) then
                        sj=j+(LTSIDXY(dr)-1)/2
                    else
                        si=i+(LTSIDXX(dr)-1)/2
                    endif
                    val=(-dtilde(LTSURFDIR(dr),1,si,sj)+dhat(LTSURFDIR(dr),1,si,sj)*SGDHAT(dr))*hcm(dr)
                    if (lo .ne. 0) call setmb2g(mtcmfd,l,k,dr,g, val,lo)
                    val=(dtilde(LTSURFDIR(dr),1,si,sj)+dhat(LTSURFDIR(dr),1,si,sj)*SGDHAT(dr))*hcm(dr)
                    call accmb2g(mtcmfd,l,k,0,indm24(g,g),val)
                enddo
            enddo
        enddo ! loop for x-direction
        ! set the real meshes in a row
        call setnxse(mtcmfd,j,nxs,nxe) 
    enddo ! loop for j-directoin
    enddo ! loop for plane
    
    call copymb2g(mtcmfdorg, mtcmfd)
    call facilu2g(mtcmfd)
end subroutine

subroutine shiftmt(ke)
    use const, only :indm24
    implicit none 
    include "hmgz.FI"
    include "xsec.FI"
    include "pspec.FI"
    real(8), intent(in) :: ke
    real(8)             :: ker
    integer             :: k,j,i,l
    
    ker = 1._8/ke
    
    call copymb2g(mtcmfd, mtcmfdorg)
    l=0
    do k=1, nz
    do j=1, ny
        do i=1, nx
            l=l+1
            if (map_exp2(i,j) .eq. -1) cycle ! for vacuum assembly
            ! update diagonal entries 
            call accmb2g(mtcmfd,l,k,0,indm24(1,1), -ker*hNuF(1,i,j)*vcm)
            call accmb2g(mtcmfd,l,k,0,indm24(1,2), -ker*hNuF(2,i,j)*vcm)
!            call accmb2g(mtcmfd,l,k,0,indm24(1,1), -ker*Xchi(1,1)*hNuF(1,i,j)*vcm)
!            call accmb2g(mtcmfd,l,k,0,indm24(2,2), -ker*Xchi(1,2)*hNuF(2,i,j)*vcm)
!            call accmb2g(mtcmfd,l,k,0,indm24(1,2), -ker*Xchi(1,1)*hNuF(2,i,j)*vcm)
!            call accmb2g(mtcmfd,l,k,0,indm24(2,1), -ker*Xchi(1,2)*hNuF(1,i,j)*vcm)
        enddo
    enddo
    enddo
endsubroutine


subroutine cmfd2gFS
!    use geom
    include "xsec.FI"
    include "hmgz.FI"
    include "pspec.FI"
    integer :: i,j,l,g
    real(8) :: lpsi
    l=0
    do j=1, ny
        do i=1, nx
            l=l+1
            if (map_exp2(i,j) .eq. -1) then
                psi(l,1)=0._8
                cycle ! for vacuum assembly
            endif
            lpsi=0
            do g=1, ng2
                lpsi=lpsi+hNuf(g,i,j)*phi(g,l,1)
            enddo
            psi(l,1) = lpsi*vcm
        enddo
    enddo
endsubroutine

subroutine cmfd2gQ(k,q)
    include "pspec.FI"
    include "arrays.FI"
    include "xsec.FI"
    real(8), intent(in) :: k
    real(8), dimension(:,:,:), intent(out):: q
    integer :: i, j, g, l, gp
    q=0
    l=0
    do j=1, ny
        do i=1, nx
            l = l+1
            if (map_exp2(i,j) .eq. -1) cycle ! for vacuum assembly
!            do g=1, ng2
!                q(g,l,1) = Xchi(1,g)*psi(l,1)/k
!            end do
            do g=1, 1
                q(g,l,1) = psi(l,1)/k
            end do
        end do 
    end do
endsubroutine

subroutine poweriter(k, iter, err, dr, ifwldt, iflog)
    use bicg2g, only : initbicg2g, solbicg2g
    use ifport
!$  use omp_lib
    implicit none

    real(8), dimension(nxy,nz)      :: psi_prv
    real(8), dimension(ng2,nxy,nz)   :: q
    real(8), intent(inout)          :: k
    real(8), intent(out)            :: dr, err
    integer, intent(out)            :: iter
    integer                         :: i, j, g, gp, maxiterin, max_iter, iti, tot_in
    real(8)                         :: errconvin, torErr, upper, upper2, lower, err_prv, ke, kd, kr
    real(8)                         :: r2, r20
    logical                         :: ifwldt, iflog
                       
!$  real(8) :: ts(3)
!$  ts(3) = 0
!$      ts(1) = omp_get_wtime()
    max_iter=iter
    iter   = 0
    err    = 1.
    dr = 1.
    if (ifwldt) kr=k
    tot_in=0                                
    ! make fission source
    call cmfd2gFS

    do while (iter .lt. maxout .and. err.gt.ccout)
        err_prv = err
        iter = iter + 1
        if (ifwldt) then
            kd=shiftk(iter, 0.01_8)
            ke=kr+kd
            k = 1._8/(1._8/kr-1._8/ke)
            call shiftmt(ke)
        endif
        ! make total source
        call cmfd2gQ(k,q)
        call initbicg2g(mtcmfd,phi,q,r20)
        if (r20<1e-10) exit
        r2=1._8;iti=0;
!        do while (r2>1e-3)  ! i=1,3
        do while (iti .lt. maxin .and. r2.gt.ccin)
            call solbicg2g(mtcmfd,r20,r2,phi)
            iti=iti+1
        enddo
        tot_in = tot_in+iti
!        
!        write(*,*) r20, r2, iti
            
        ! make fission source from updated flux
        psi_prv = psi
        call cmfd2gFS
                
        ! evaluate k                                            
        upper = 0.
        lower = 0.
        upper2= 0.
!!$omp parallel do private(i), reduction(+:upper,lower,upper2)       
        do i=1, nxy                                                
            upper = upper + psi(i,1)*psi(i,1)
            lower = lower + psi(i,1)*psi_prv(i,1)
            upper2 = upper2 + (psi(i,1)- psi_prv(i,1))*(psi(i,1)- psi_prv(i,1))
        end do
!!$omp end parallel do
        if(ifwldt) then
            k=k*upper/lower
            kr=(k*ke)/(k+ke)
        else
            k=k*upper/lower
        endif        
                
        ! evaluate err 
        lower = sqrt(upper)
        upper = sqrt(upper2)
        err = upper/lower
    end do
!$      ts(2) = omp_get_wtime()
    if (ifwldt) k=kr
    dr = err/err_prv  
    if(iflog) then
        write(*,*) iter, err, tot_in
        write(*, '(" Solution from FDM : k=",f7.5,1x,"dr=",f7.5)'), k, dr
        write(*, '(" Computing time    : ", f12.3,1x,"(sec)")') ts(2)-ts(1)
    endif
    
    contains
    
function shiftk(iter, delk) result (ret)
!    real(8) :: SHIFT_GUIDE(17)
!    data SHIFT_GUIDE /1.0_8, 0.75_8, 0.50_8, 0.35_8, 0.25_8, &
!                    0.20_8, 0.15_8, 0.125_8, 0.100_8, 0.075_8, 0.050_8, &
!                    0.035_8, 0.025_8, 0.020_8, 0.015_8, 0.0125_8, 0.0100_8/
    real(8) :: SHIFT_GUIDE(7)
    data SHIFT_GUIDE /1.0_8, 0.50_8, 0.25_8, 0.125_8, 0.075_8, 0.035_8, 0.01_8/
    integer, intent(in) :: iter
    real(8), intent(in) :: delk
    real(8) :: ret
    ret = delk
    if (iter .le. 7) then
        if (ret<SHIFT_GUIDE(iter)) ret = SHIFT_GUIDE(iter) 
    endif
end function
    
end subroutine


endmodule
